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Abstract. Laplace transforms which admit a holomorphic extension to some sector strictly 
containing the right half plane and exhibiting a potential behavior are considered. A spectral order, 
parallelizable method for their numerical inversion is proposed. The method takes into account the 
available information about the errors arising in the evaluations. Several numerical illustrations are 
provided. 
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1. Introduction. In a variety of situations, the problem arises of inverting nu- 
merically the Laplace transform U(z) of a given mapping of interest u(t). Roughly 
speaking, it turns out that the wider the set W where U(z) can be computed is, the 
easier the inversion results. For instance, if W is an interval (a, b) then the numerical 
inversion becomes an ill-posed problem ^IHlEj. On the other hand, if W is the 
complement of some bounded region, then the efficient Talbot's method I25j is at 
hand. 

In the present paper we focus on the particular situation where W is a sector 
symmetric with respect to the real axis, strictly containing the right half plane, and 
we assume that U(z) exhibits a potential behavior on W. We say then that U(z) 
is sectorial. Precisely, there is a renewed interest in the numerical inversion of sec- 
torial mappings |1U1 111! 1141 1161 |2"T] , mainly due to its applicability to linear, non- 
homogeneous evolution equations of parabolic type (both in the context of abstract 
IVP's and Volterra equations), as well as their discretizations in space [21 EH- Notice 
that the applicability of the inversion approach, in the sectorial setting, demands in 
practice that the source term of the parabolic equation must be approximated effi- 
ciently (at least locally) by holomorphic mappings This difficulty is overcome 
in |121 12U| , where the ideas in the present paper are adapted so as to provide accu- 
rate reconstructions of the traditional Runge-Kutta approximations to the solutions 
of such parabolic problems. These reconstructions require no regularity on the source 
term of the problem. 

In the present paper we consider the issue of the numerical inversion of sectorial 
mappings by itself. To fix ideas, let 



be a locally integrable mapping, taking values in a Banach space X, with exponential 
growth. Denote by 
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u : (0,+co) -> X 




its Laplace transform. We will always assume that U(z) admits a holomorphic exten- 
sion to the complement W of some acute sector 

E 4 = {zeC: \axg(-z)\ <8}, < 5 < f > 
an that there exist constants M > and \i G M such that 

\\U{z)\\<^, zi^s- (1-2) 

The last requirement, with fi > 1, means that u admits a bounded and holomorphic 
extension to any sector of the form | arg(z) < 8', with < 8' < tt/2 — 8. If \i < 1, we 
select an integer number m > 1, with m + fi > 1, and set = f7 (z)/z m . Then, by 
the previous remark, V(z) is the Laplace transform of a mapping t> : (0, +oo) — * X, 
which admits a bounded and holomorphic extension to sectors with semi-angle 8' as 
before, and now u is understood to be the derivative of order m of v. 
Notice that in case U (z) satisfies a similar inequality 

M 

\\U{z)\\ < z£oj + Z s , 

for some lu G R, then, by using the shifting theorem, the inversion of U(z) is reduced 
to the one of a Laplace transform U(z) fulfilling i|1.2[l . Since the respective originals 
u{t) and u(t) are related by u(t) = e ult u{t) 1 then we can just approximate u(t). This 
is why the analysis is restricted to the situation to = 0, i.e. to (|1.2fl . 

The goal is to numerically reconstruct u from knowledge of a moderate number of 
evaluations of U(z) at suitable nodes z £ £5. Let us point out that, from a practical 
point of view, it is essential to take into account that these evaluations are going to 
be affected by errors. 

The starting point of the method we propose is the well-known inversion formula 



u(t) = — [ e tz U(z)dz, t>0, (1.3) 

where L is a suitable path connecting —zoo to +ioo which, in our setting, can be 
taken so as to guarantee the absolute convergence of the integral appearing in 1|1.3|) . 
As in |1(J1 1141 IKij , we choose L the branch of a hyperbola parametrized by a mapping 
S : (—00, +00) — > C admitting a holomorphic extension to a horizontal strip around 
the real axis. The numerical method we propose is simply the truncated trapezoidal 
rule, applied to the definite integral arising after parametrizing by S, used with 
2n + l nodes Xk = kh, —n < k < n, and a suitable step size h > 0. The properties of S 
allow us to use the ideas and results in |231 , where the trapezoidal rule applied to 
holomorphic mappings on strips is considered. Let us comment that the fast decay of 
our integrand yields an improvement of the more general estimates in |2,'ill2 lj . 

Very often, for instance in the context of IVP's (see Illustration |3 in Section |3J), 
the main computational effort of the method is due to the evaluations of U(z) at the 
nodes Zk — S(xk), —n < k < n. An important feature of the present approach is 
that the same evaluations can be used to approximate u(t) at different t > |14l 
Accordingly, our goal is to obtain a uniform error estimate for the approximation 
of u(t) on intervals of the form [to, Ato], with given to > and A > 1, rather than 
at a fixed t > 0. Essentially, this was the aim in whose basic estimates we 
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borrow. Notice also that the algorithm presents two levels of parallelism since, first, 
the evaluations of U (z) at the involved nodes and, second, the evaluations of u(t) at a 
selected finite set of values of t £ [to,Ato], can be carried out on different processors. 

In the present paper, by considering a different choice of the geometrical and scale 
parameters from the one in |14) . we improve the results there in two different ways: 

(i) We get a better error bound, which now turns out to be a genuine spectral 
estimate of the form 0(e~ cn ). 

(ii) We also get a weaker dependence of the exponential factor c on A, since now 
c = (9(l/lnA). 

This means, in practice, that with a moderate number of evaluations of U (z) we can 
accurately approximate u[t), uniformly on intervals [to,ti] with A = t±/to » 1, let 
us say A = 50. 

On the other hand, for the choice of parameters we propose, the precision p used in 
the evaluations of U(z) at the required nodes plays a more relevant role than in |14| . In 
fact, ignoring that we always have p > would result in large actual errors for n >> 1, 
as simple numerical experiments show (see Illustration I in Section^. This drawback 
is overcome by minimizing the estimate we get for the actual error (Theorem^, which 
leads to a (p, n)-dependent choice of parameters. With this choice, the actual error 
finally behaves for moderate n like 0(e~ cn ), with c = 0(1/ In A), and for large n like 
O(p). This optimal choice of parameters demands, of course, some information about 
the size of p. In the absence of it, we propose an n-dependent choice of parameters for 
which the actual error behaves like 0(p + e~ cn ), with c = 0(l/(lnn + In A)). All the 
above estimates are uniform on to < t < Ato, with fixed to > and A > 1. Moreover, 
the error constants are made explicit in the analysis and turn out to be reasonable. 

The outline of the paper is as follows. In Section |3 we describe the numerical 
method and show, in Theorem^ how to achieve (i) and (ii). The propagation of 
errors is studied in Sectional The choice of parameters is considered in SectionQ]and 
four simple numerical illustrations of the theoretical results are provided in Sectional 

2. The numerical method. Given 5 in Ijl.f I) and following the ideas in |14|. 
we select a, d > such that 

0<a-d<a + d<^-5 . (2.1) 

Defining 

T(w) = 1 - sin(a + iw) (2.2) 

this mapping transforms each horizontal straight line Imw = y, —d < y < d, into the 
left branch of the hyperbola given by 



Re z — 1 \ 2 / Im z 



sin(a — y) J \cos(a — y) 



= 1, (2.3) 



with center at (1, 0), foci at (0, 0) and (2, 0), whose asymptotes make angles ±[7r/2 — 
(a — y)] with the real axis. Therefore, T transforms the horizontal strip 

D d = {z e C : |Imz| < d} 

into the region in the complex plane limited by the left branches corresponding to 
y = ±d in l|2~3jl. 



Introducing a parameter A > 0, the parametrization of T in can be denned 

as 

T = {XT(x) : x e K}, 

i.e. r is the branch of a hyperbola corresponding to the image of the real axis under 
S = XT. This results in 

/+oo 
G t (x)dx , t>0, 
-oo 

where Gt '■ Dd — * X, t > 0, is the mapping 

G t (w) = -^exp(\tT(w))U(\T(w))T'(w). 

2lTl 

Once the parameters a, d, and A have been fixed, we set Xk = kh, k € Z, and consider 
the approximation to u(t) given by 

n 

u n (t) = hJ2 G t(xk), t > 0. (2.4) 

k— — n 

The proof of the main result in ^5] (Theorem 2), shows that for \i = 1 in (| 1.2(1 

\\u(t)-u n (t)\\<MM^d)^tsin(a-d)).e^[ e2 J^ (2.5) 
where 



, „ 2 /l + sin(a + rf) 

<p(a,d) = -J- — — — , 

7r y 1 — sm[a + d) 

and L(x), x > 0, is the function 

L(x) = l+|ln(l-e- a: )|. 

Notice that L(x) is decreasing in x, L(x] ~ | In ar | , for x — » + and £(a;) tends to 1, 
for x — > +cxd. 

As we commented in the Introduction, in many applications the computational 
effort to obtain u n (t) is mainly due to the evaluations of U(z) at z = AT(xfe), — n < 
k < n, but these evaluations could be carried out in parallel. Another attractive 
feature of (|2.4Jl is that the same evaluations of U (z) can be used to compute u n (t) for 
different t > 0. In fact, as we see below, with the appropriate choice of parameters, 
we can use the same evaluations of U (z) so as to have a spectral estimate 

\\u{t)-u n {t)\\ =0(e~ cn ), 

uniform on intervals to <t<t\. The exponential factor c turns out to depend weakly 
on the ratio A = ti/to, given that c = 0(l/lnA). 

For simplicity the next theorem is restricted to the situation /x = 1 in l|1.2(l . The 
cases fi > 1 and /i < 1 are treated in subsequent remarks. 

Theorem 1. Assume that U satisfies Hl.fy) with fj, = 1. Fixing a and d according 
to \2.1[) . for to > 0, A > 1, < 9 < 1 cm<i n > 1, the following choice of parameters 
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with 

a<9) = arccoshf- -r- . , . 

\(1 — v) sin a/ 

yields the uniform estimate on to < t < A<o 

\\u(t) - u n (t)\\ < M ■ <p(a, d) ■ L(\t sm(a - d)) ■ , (2.7) 

where 

2nd 



e n (B) = exp ( 



a(6) 



n 



The theorem shows, just by selecting any < 9 < 1, a genuine spectral order of 
convergence in n of the form 0(e _CTl ), where c = (9(l/lnA) (cf. jl()l I14| ). 
Proof. Set a = Xt . For t < t < At , 12. ty implies the uniform bound 

IK*) - u n (t)\\ < M ■ <p(a, d) ■ Lip sin(a - d)) ■ _ 1 + egrina Lh(n/>) ) ' 

Our choice of h and A is precisely the one guaranteeing that 
cxp ( — — ^ = exp(cr sin a cosh(nh)) 



h J v " e n (ey 

hence 

2e -2*d/h 2e n {6) 



< 



g2Trd/h ^ g(T sin a cosh(n/i) — ^ ^ — 2ird/h ^ Cn(^) 

The proof ends after remarking that 

e A °e n (0) = e n {9) e - l e n {6) = e n (9) e . □ 

To end the section we comment, in the two following remarks, on the situation 
fi ^ 1 in 11. 2f) . We omit details in the proofs, which are completely analogous to the 
one of Theorem ^ 

Remark 1. Assume that U satisfies (|1.2|) with fx > 1. By Remark 1 in we 
have 

e At / 1 1 \ 

IK*) - Un(t) || < M • <pia, d, /X) • L(At Sin(a - d)) • ^ZT { e27ld/h _ 1 + e A t sinaco S h(n fe ) J ' 

where 



, , , 2 l + sinfa + d) 
(p(a, d, n) - 



7r V (1 — sin(a + d)) 2 ^ 1 ' 

Thus, for < 9 < 1, the same choice of values for h and A as in Theorem ^ gives the 
bound 

IK*) - «n(*)|| < M ■ <pia, d, fi) • L(Xt sm(a - d)) ■ A 1 "" ■ ^"^1 , 
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uniformly for i < t < Ai . This estimate is again spectral in n, since 

Remark 2. Assume now that U satisfies (|1.2|) with fx < 1. By Remark 1 in |T3|. 
for a fixed s G (0, 1), there holds 



e xt ( 1 

\\u(t)-u n (t)\\ < M-(p s (a,d, ^■L(s\tsm{a-d))-- r -j I ^ 27rd/h _ - 
where now 



, , . 2 / 1 + sin(a + d) ( N 1 "** 



7r y 1 — sin(a + d) \ (1 — s)e sin(a — d) 
In this situation, for 9 £ (0, 1) we choose 

1 2irdn(l - 9) 

h=-a s {9), x =—r-\ — 7^> 

where 

a s (9) = arccoshf— ^— 

Vs(l - 9) sin a 

Setting 

(-2-Kdn 

£ s ,n{v) = exp 

we get the spectral estimate 



a a (0) /' 



2e (9)^ 

\\u(t) - u n (t)\\ < M ■ (p s (a,d,fJ,) ■ L(s\t sm(a - d)) ■ t%~ 1 - — — — j-rr, 

uniformly for to < t < Aio- 

3. Error propagation. Numerical experiments (see Section [SJ), show that for 
large values of n the estimate (|2.7|l is not longer true in practice. The explanation 
of this apparently contradictory behavior lays in the influence of the errors when 
evaluating U and the elementary functions involved. For the sake of simplicity, we 
consider first the case fi = 1 in (|1.2|l . The situations /i > 1 and // < 1 are considered 
in subsequent remarks. 

Let Zk = \T(xk), —n < k < n, be the nodes used in (|2.4|) . Clearly, in practice, 
as numerical approximation to u(t) we actually obtain 

n 

u n {t) = u k(f)U k , (3.1) 

k— — n 

where, for — n < k < n, tOk{t) 6 C and Uk £ X, are approximations to 

-7T- exp(Atz jt )T'(a; fc ), 
2m 



and U(z k ), respectively. 

To estimate the actual error \\u(t) — u n (t)\\ we need to make some assumptions 
on the approximations used. To this end, we are going to focus on two frequent 
possibilities, depending on whether we have information on absolute or relative errors 
due to the evaluations. To be precise, we are going to assume that there exists p > 
such that, simultaneously for all —n < k < n, we have either 

\\U{z k ) - [/ fe || < p and w fc (t) = exp(Xtz k )T'(x k ) (3.2) 

2tti 



|| exp(Xtz k )T'(x k )U(z k ) - u k (t)U k \\ < p|| exp(\tz k )T'(x k )U(z k )\\. (3.3) 

Situation (|3.2|l arises for instance when U k ~ U(z k ) are provided by means of 
some auxiliary routine, let us say by solving a linear system, with prescribed accuracy 
p and moreover the errors due to the evaluations of the elementary functions involved 
turn out to be negligible compared to p. Situation l|3.3|) is typical when U(z) is an 
elementary function. 

The next theorem yields an estimate of the actual error for these situations. We 
maintain the notation introduced in Theorem ^ 

Theorem 2. Assume that U satisfies with fi = 1. Fix a, d according to 

\2.1\) . For to > 0, A > 1, < 9 < 1 and n > 1, select the 'parameters 

h=-a(0), A = ■ 

Assume also that w k (t) G C, to < t < t\, U k £ X, — n < k < n, satisfy either 13.2(1 
or Then, the actual error is estimated by 

IK*) - u n (t)|| < M ■ $ • Q ■ [se^ef- 1 + ^2(0) ) ' (3 - 4) 

uniformly on to < t < Ato, where either 
(a) e = p/(Mt ), 



.2 /l + sin(a + d) 1 
$ = max < - < ' 



7r V 1 — sin(a + d) ' ire sin a 



and 



in case \3. 2\) holds, or 
(b) e = p, 



{In n 
2L(Xta sm(a — a)), - 
In n — 1 



In n ( A^o sin a 
L 



2n V In n 




1 + sin(a + d) 
tt y 1 — sin(a + d) 

and 

Q = max{2i(At sin(a - d)), l/2(h + L(\t sin a))}, 



in case J^3.y\) holds. 
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Notice that Q depends logarithmically on a, d, 1 — 8 and A. 

The estimate (|3.4f) given by the theorem, with a fixed < 8 < 1, shows again a 
spectral order of convergence 0(e~ cn ), with c = 0(1/ In A), but only for moderate n, 
to be more precise, as long as e„(6>) > e. In fact, for fixed 8, l|3.4|) goes to +oo as 
n — > +oo. However, this apparent drawback is overcome by selecting 8 in a suitable 
way, as we explain in Section 4. 

Proof. By writing 

\\u(t) - u n (t)\\ < \\u(t) - u n (t)\\ + \\u n (t) - u n (t)l 
and noticing that, for the corresponding Q, 1|2.7|1 implies 

\\u{t)-u n {t)\\<M^-Q en{6)(> 

1 - e n (B) 

the proof is reduced to show that 

K(t) - u»(*)|| < M • $ • Qee n {8) e -\ (3.5) 

Assume first that Il3.2[l holds. This situation was already studied in JJ], where 
it was proved that 



||«n(t)-«n(t)|| < 



plan 



27re(lnn — 1) sin a t 



In n 



2L 



Xt sin a 
Inn 



whence, after recalling that e = p/(toM) and noticing that 



„AAtn 



£n(0) 6 



(3.6) 



(3.7) 



we readily obtain l|3.5l) . 

Assume now that (|3.3|) holds. Proceeding as in the proof of Lemma 1 and Theo- 
rem 2 in and denoting 



( t\\ 2 /1 + sina 
7r V 1 — sin a 



we get 



IMt) -u„(*)|| < 



pMe 
2tt 



< 



At ra 

E 

k=—71 

M<p{a,Q) w 



-Ais 



k cosh : 



T'(x fc ) 



T(x k ) 



pe xt h z 



Xt sin a cosh a:^ 



fe=-T 



My(a,Q) At 

— o " 



-Xt sin a cosh x 



dx 



<M^O) pe Xt {h + L{Xtsina) y 

Hence, using again 13.7|l and the inequality (p(a, 0) < (f(a, d), we deduce (13.511 . □ 
The behavior p ^ 1 in (|1.2I) is considered in the following remarks, whose proofs 

are a combination of Remark 1, Remark 2 and the arguments used in the proof of 

Theorem in |14| . Notice that (|3.6|l is independent of p. 

Remark 3. Assume that p > 1 in (|1.2|l and fix < 8 < 1. Then, for the choice 

of parameters in Theorem |21 and uniformly on to < t < Atg, we have: 
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(a) in case l|3.2|) holds 

\\u(t) - u n (t)\\ <M-$-Q- ( :< n U-))" ' A 

with e = p/(Mt ), 



1 , e n(Q) 



1 - e n {9) 



$ = max < — 



1 + sin(a + d) 



7r V (1 — sin(a + d)) 2 ^ 1 ' ne sin a 



and 



Q = max ^2L(A^o sin(a — d)), 



In i 



In n — 1 



In n ( Xto sin a 



2n 



Inn 



(b) in case <|3.3|) holds 

IK*) ~ <M$Q- A 1-M ■ (ee^e)"- 1 + ^"^(g) )' 



with £ = p, 



and 



2 / 1 + sin(a + d) 
(1 - sin(a + d)) 2 ^" 1 



Q = max{2L(At sin(a - d)), l/2(/i + L(At sin a))}. 

Remark 4. Assume that u < 1 in II. 2[) and fix < s, < 1. Then, for the choice 
of parameters in Remark [21 and uniformly on to < t < Aio, we have: 
(a) in case (|3.2() holds 



K*) - «n(*) II < M ■ $ • Q • £e S) „(^) ( '- 1 + fg 



with s = p/(Mt ), 



,2 1 + sin(a + d) 
<£ = max < - 1 



1-fi 



-1 £s,n{9) B 
l-e a ,nW>/' 



1-M 1 



7r y 1 — sin(a + d) V (1 — s)e sin(a — d) / ' 7re sin a 



and 



Q = max < 2L(s\to sin(a — d)), 



In i 



In n — 1 



In n / Aio sin a N 
2n \ Inn / 



(b) in case l|3.3l) holds 



IK*) - un(*)|| < M • $ • Q • ^-"ee,,,,^) - 1 + tff -1 - 



e s ,n(^) e 
-e»,n(e) 



with e = p, 



and 



$ = 



2 /l + sin(a + d)/ 1 - /x x 1 -^ 



7r y 1 — sin(a + d) V (1 — s)e sin(a — d) 



max{2L(sAio sin(a — d)), l/2(h + L(sXto sin a))}. 
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4. The choice of parameters. With Theorem|2in mind, we now try to derive 
a strategy for the choice of parameters. First of all, (|3.4|) shows that it is of interest 
to select a away from zero and a + d away from tt/2. The dependence of the actual 
error on a — d is less important, since it is logarithmic. 

Suppose a and d have been already chosen, then for a given n we take h and A as 
indicated in Theorem|21and we fix < 9 < 1. Assume also that we have an estimation 
of p and set e = pj (Alio) or e = p as in Theorem 2. Then, since in practice we always 
have p > and hence e > 0, it turns out that ee n (9) 6 ^ 1 — > +oo, as n — » +oo. Hence, 
it is clear that increasing the number of nodes might result in a worse estimate l|3.4[) . 
In fact, increasing n may result in worse approximations, as Illustration^^ Sectional 
shows. 

To overcome this drawback we let 9 be a free parameter for the moment. Given 
e > and n, after selecting a and d, neglecting the logarithmic factor Q and taking 
into account that typically e n (8) << 1, the best thing we can do is to choose < 9 < 1 
so as to minimize the term 

senie) 6 - 1 + e n {9) 6 (4.1) 

i.e. we must tune 9 depending on e > and n. By a direct calculation it can be 
proven that the first derivative of e n (9) with respect to 9 is increasing in 9. The 
same is true for e n (9) e (in this case the proof, though elementary, is more difficult). 
We conclude that the expression in (|4.1(l is a convex function of 9. Moreover, its limit 
either for 9 — » 0+ or 9 — > 1— is +oo. Therefore, <|4.1|) attains its minimum exactly for 
one value d E<n G (0, 1), which is the one we propose to be used. Though it is not easy 
to express the dependence of 9 £ n on n and e, this can be easily done numerically (see 
Section 5). 

Since, up to logarithmic factors, the choice 9 — 6 £>n in .4|1 is optimal, it is clear 
that with this choice we get for the actual error: 

(a) A spectral order of convergence 0(e~ cn ) with c = 0(l/lnA), for moderate 
values of n, since this is true for any value < 9 < 1. 

(b) The errors are not propagated. In fact, already with the non-optimal choice 



reads 

\\u(t) - u n (t)\\ -0(e + e- c "), (4.2) 

uniformly on to < t < Aio, with c = 0(l/(ln A + Inn)). This remark tells us 
that, for large values of n, the actual error saturates at level e, as observed 
in the numerical experiments (see Section 5). 
In the previous discussion it was essential to assume that we had some information 
about e. Notice that, even in case we do not have such an information, the choice 
9=1 — 1/n, which led to 1)4.2)) . is always available. This bound is almost spectral in 
n, depends weakly on A and prevents error amplification. 

5. Numerical illustrations. In this section we give four numerical illustrations. 
The first two ones concern elementary Laplace transforms which are assumed to be 
computed with a relative error of order p w eps, where eps stands for the machine 
precision (eps = 10 -16 in our computations). In the last two illustrations we do not 
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Number of nodes n 



Fig. 5.1. lnmax t g[ t0iJ \t ] ||"(*) — u n(t)\\ versus n for u in Illustration QJ with 8 = 0.5 fixed, 
a = 0.7, d = 0.6 and to = 1. The gray line corresponds to A = 50 and the black one to A = 5. 

assume any information about the errors due to the computations of the Laplace 
transforms. 

Illustration 1. We first show by means of a simple example, that for n >> 1 
<|2.7|l fails in the presence of errors in the evaluations. To this end, we consider the 
mapping u(t) = e - *, whose Laplace transform is U(z) = 1/(1 + z). 

This function satisfies ljl.2jl for all 6 > and M = l/sind. We fix 9 = 0.5, a — 
0.7, d — 0.6 and choose the parameters h, X as stated in the theorem for all the values 
of n. In Fig. I5.1l we plot in a semilogarithmic scale the absolute actual error, i.e. 

In max \\u(t) — u n (t)\\ 
te[t ,A*o] 

versus n (recall that u n (t) stands for the actual computed approximation to u(t), 
see P-lf) ). This is done for A = 5, 50 and to = 1. This figure shows that the error 
decays exponentially for the first values of n, saturates near e level and then grows 
like 0(e cn ). 

Next we tune parameters as explained in Section 4. For A = 5, 50, in Fig. 15.21 
(left) we plot the optimal values of 9 against n. In Fig. 15.21 (right) we plot 

In max \\u(t) — u n (t) || 
te[t ,At ] 

(continuous line) and the logarithm of the corresponding values of the theoretical 
error estimate (dashed line) obtained in Theorem [3 versus n, once 9 is optimal. We 
maintain a — 0.7, d — 0.6 and to = 1. 

Illustration 2. Take [3 = 1.5 and set 
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Fig. 5.2. Left: Optimal 8 versus n. Right: Natural logarithms of max t g[ tn J \t n ] ||tt(t) — u n (£)|| 
(continuous) and the theoretical estimate (dashed) versus n, for u in Illustration^! The gray lines 
correspond to A = 50 and the black ones to A = 5. 

i.e., U(z) is the Laplace transform of 

u(t) = Mpi-t?), 

where Mp stands for the Mittag-Leffler function of order j3 (see JE])- Notice that 
U satisfies (Ol) for any 8 E (w/3, n/2), with ^ = 1 and M = l/sin(/3(7r - 5)). 
We consider here as exact solution the one computed with 500 nodes and take a = 
7r/12, d = 0.25 and t = 1. 

This example was already considered in |14j . In order to compare the performance 
of the strategy proposed in with the one proposed in the present paper, we first 
compute u n {t) by selecting the parameters as in I n Fig- E5H (left) we plot in 

semilogarithmic scale the theoretical estimate and actual errors for A = 2, 5, which 
are acceptable. In Fig. 15.31 (right) we do the same for A = 50 and conclude that the 
approach in ^3] is n °t at all useful for large values of A. However, the corresponding 
computation by using the strategy in Section 4, yields the plot in Fig. 15.41 which 
shows a satisfactory spectral order of convergence even for A = 50. 

Illustration 3. We consider the inhomogenous heat equation on the unit square 
Q = (0, l) 2 with zero initial value and a convective heat flux at the boundary 

u t (t,x) = Au{t,x) + f(x), for x e Q, t > 0, 
d u u{t,x) = —u(t,x), for x G dil, t > 0, (5.1) 
u(0, x) = 0, for x e O, 

where / is the indicator function of the rectangle R — [0.6, 0.8] x [0.2, 0.8], i.e. / = 1 
on R and / = elsewhere. 

Problem l|5.1|) is semi-discretized in space by using linear finite elements on a 
triangular grid. Denoting by Vh C L 2 (tt) the space of elements and by Uh(z) the 
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Fig. 5.3. Natural logarithms of max t gr tch Atn] ll M ffl — II (continuous) and the theoretical 

estimate (dashed) versus n, for u in Illustration lH vroceedina as in [TTJ for 8 = tt/3, t = 1. The 
gray lines correspond to A = 50 and the black ones to A = 5. 



Laplace transform of the semi-discrete solution Uh{t), we get 

U h (z) = -(z-Ah^Phf, 

z 

with A/j : Vh — » V/, the discrete Laplacian and Ph the orthogonal projection of / onto 
Vh- Now, for fixed h > 0, we try to approximate it^t) by inverting Uh(z). Notice that, 
since A/j is definite negative, certainly Uh{z) satisfies (ll.2|) for any < S < tt/2 and 
M = l/(/x/ 1 sin(5)), with —Hh the highest eigenvalue of A/j. Notice also that, working 
in coordinates relative to the standard basis of elements, Uh(z) is represented by a 
vector valued mapping \Jh(z) satisfying 

zM h XJ h (z) + S h V h (z) - -f h , 

z 

where Mh and Sh stand for the mass and stiffness matrices and where f/j is the vector 
formed by the scalar products of / with the elements of the basis. Thus, one evaluation 
of U(zk) at a given node z^ requires the solution of one linear system of the above 
form. 

In the experiment we generate a mesh, shown in the left of Fig 15.51 with 542 
triangles by means of the mesh generator Triangle |22j . Linear systems are solved 
using MATLABs sparse LU factorization UMFPACK. Since Uh(t) is unknown, the 
errors are estimated in the £ 2 (f2)-norm with respect to a reference solution w/^soo^) 
obtained with 500 nodes. In the absence of precise information about p, both for 
this reference solution and for the rest of the approximations Uh,n{t) to Uh(t), we tune 
9 = 1 — 1/n, as indicated in Section^] In Fig. 15. 61 for the parameters a = 0.7, d = 0.6, 
t = 0.01 and = 1 - 1/n, we plot lnmax te|t(ljAfo ] ||%,50o(<) - u h , n {t)\\ against n, for 
A = 5, 50. This plot shows the predicted behavior. 
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Fig. 5.4. Natural logarithm of max f £\t n ,At n ] \\ u (t) ~ u n (t) \\ (continuous) and the theoretical 
estimate (dashed) versus n, for u in Illustration^^ The gray lines correspond to A = 50 and the 
black ones to A = 5. 
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Fig. 5.5. Left: Mesh ofQ, with the set Ft indicated by dark-gray. Right: Temperature distri- 
bution at t = 0.5 in false-color representation, (white corresponds to temperature 1 and black to 
0) 

Illustration 4. We consider again the Laplace transform U(z) = 1/(1 + z) of 
the exponential function u(t) — e~* as in Illustration^ The values of a, d and to are 
again 0.7, 0.6 and 1, respectively. 

We add on purpose perturbations of maximum size 10 -4 to the evaluations of U 
at the required nodes. Thus, we use (|3 . 1ft with 

U k = U(z k ) + rj k , -n < k < n, 

with \rjk\ < p = 10~ 4 . Now we try to approximate u(t) without using the available 
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Fig. 5.6. Left: Natural logarithm o/max t g[^Ato] II^W - (continuous) and the theoretical 
estimate (dashed) versus n, for u in Illustration]^ The gray lines correspond to A = 5 and the black 
ones to A = 50. 



information about p. In this situation, as explained in Section^ we take 6 = 1 — 1/n. 
In fact, we compare two types of perturbations: 

We first generate complex, random, independent perturbations in such a way 
that \r]k\ and arg(^fc) are uniformly distributed on [0, 10 -4 ] and [0, 2ir], respectively. 
In Fig. 15.71 (left), we show the resulting actual error, which behaves much better than 
predicted by l|4.2(l . The explanation is that cancellations are likely compensating the 
effects of the independent random perturbations. A finer analysis of the observed 
behavior is out of the scope of the present paper. 

Secondly, for each —n < k < n, we consider the perturbation 

?? fc = ICr 4 cxp(-iarg(cj fc (t ))), 

with uJk{to) defined in [|3.1JI . These perturbations correspond to the worst possible 
case in (|3.2f) . for t = to = 1. Now, the resulting actual error, plotted in Fig. 15.71 
(right), fits quite well with (|4.2|) . 
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